close all
clear
clc

% Common value for simulation
c0 = 299792458; % light speed
nmax = 1; % free space
LAMBDA = 0.08; % wavelength
freq = c0/LAMBDA;

Zlength = 2.5; % total Z for simulation

% =================================================================
% Start calculate grid resolution
NRES = 10; % resolve the wave with at least 10 cells, better >=10
NDRES = 4; % normally 1~4 for resolution of feature size
dmin = 0.1;

 % Compute default grid resolution
dz1 = min(LAMBDA)/nmax/NRES; 
dz2 = dmin/NDRES; % dmin means the smallest feature size
dz = min(dz1,dz2);
% Snap grid to critial dimensions
N = ceil(Zlength/dz);
dz = Zlength/N;

Nz = floor(Zlength/dz);
% End calculate grid resolution
% =================================================================

% =================================================================
% Start calculate delta_t and tau 
Nt = 10; % Nt >=10
tau = 0.5/freq; % tau <= 1/(pi*freq) 

% Delta T resolution
dt1 = tau/Nt;
dt2 = dz/(2*c0);
% if using absorbing boundary, need make sure dt = dz/(2*c0),
% in order to make sure each cell need two time step

dt = min(dt1,dt2);
% End calculate delta_t and tau 
% =================================================================

% =================================================================
% Start calculate simulation time
t0 = 6* tau; % at least 6 tau for a whole pulse
tprop = Nz*dz*nmax/c0;
T1 = 12 * tau; % total simulation time should >12 tau
T2 = 5 * tprop; % total simulation time should > 5 bounces
Ttotal = T1 + T2; % for highly resonant devices, this is not enough
% End calculate simulation time
STEPS = round(Ttotal/dt);
% =================================================================


% Initialize Materials to free space
ER = ones(1,Nz);
UR = ones(1,Nz);

% Compute updated coefficients
mEy = (c0*dt)./ER/dz;
mHx = (c0*dt)./UR/dz;

% Initialize Ey and Hx to zero
Ey = zeros(1,Nz); 
Hx = zeros(1,Nz);


E3=0; E2=0; E1=0;
H3=0; H2=0; H1=0;

fig=figure;
set(fig,'Name', 'FDTD 1D Simulation');
set(fig,'NumberTitle', 'off');

% Main FDTD Loop
for T = 1 : STEPS
    
    % Update H from E
   for nz = 1 : Nz-1
       Hx(nz) = Hx(nz) + mHx(nz)*(Ey(nz+1) - Ey(nz));
   end
   H3=H2; H2=H1; H1=Hx(1);
   % Dirichlet Boundary Conditions
   % Hx(Nz) = Hx(Nz) + mHx(Nz)*(0 - Ey(Nz)); 
   % Absorting Boundary Conditions, pay attention to delta time condition
   Hx(Nz) = Hx(Nz) + mHx(Nz)*(E3 - Ey(Nz)); 
   
   % Dirichlet Boundary Conditions
   % Ey(1) = Ey(1) + mEy(1)*(Hx(1) - 0); 
   % Absorting Boundary Conditions, pay attention to delta time condition
   Ey(1) = Ey(1) + mEy(1)*(Hx(1) - H3); 
   for nz = 2 : Nz
       Ey(nz) = Ey(nz) + mEy(nz)*(Hx(nz) - Hx(nz-1));
   end
   E3=E2; E2=E1; E1=Ey(Nz);
   
   % Inject the source
   t = T * dt;
   g=exp(-((t-t0)/tau)^2);
   Ey(30) = Ey(30)+g;
   
   %Plot
   plot(Ey,'-b','LineWidth',2);
   hold on
   plot(Hx,'-r','LineWidth',2);

   title(sprintf('Step: %d',T));
   xlim([0 Nz])
   ylim([-1.5 1.5])
   drawnow;
   hold off
   pause(0.01);
end



